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Preface 


The  purpose  of  this  study  was  to  investigate  the  application  of 
algorithms  derived  from  the  study  of  fractal  geometry  to  the  general 
problem  of  pattern  recognition.  The  original  goal  was  to  decompose  an 
arbitrary  input  scene  into  a  fractal  basis  set.  It  was  hoped  that  this 
basis  set  would  be  orthogonal,  such  that  expansions  similar  to  Fourier 
series  oould  be  derived.  The  problem  proved  more  difficult  than 
anticipated,  and  the  application  of  fractal  geometry  to  image 
segmentation  was  then  investigated  as  groundwork  for  further  study  on 
the  basis  set  problem.  Two  widely  used  methods  for  calculating  frcctal 
dimension  and  an  original  hybrid  technique  were  investigated.  All 
three  techniques  appeared  to  be  effective  to  3ome  extent  in  segmenting 
Images,  but  shared  a  common  problem  with  the  establishment  of  suitable 
thresholds  for  robust  segmentation. 

Through  the  course  of  this  reserved  time  of  study,  I  have  had  a 
great  deal  of  encouragement  and  assistance  from  those  around  me.  I 
would  like  to  thank  my  faculty  advisor,  Dr.  H.  Kabrisky,  for  his 
receptiveness  to  this  topic  of  study,  Haj.  P.  Ambum  for  his  enthusiasm 
and  energy,  Capt.  Steve  Rogers  for  his  constructive  comments  as  a 
committee  member  and  reader,  and  Capts.  Richard  Roberts  and  Dennis  Ruck 
for  their  assistance  with  the  ITEX^TOJiMS  conversion  program.  Ky 
family  has  been  a  constant  source  of  relaxation  and  encouragement,  and 
t  wish  to  thank  them  for  their  understanding  when  studies  beckoned,  and 
their  prayers  for  Godspeed  and  success.  I  am  particularly  indebted  to 
my  lovely  wife  Jennifer  tor  her  unending  encouragement  through  those 
times  when  answers  seemed  distant;  and  my  precious  daughter  Emily,  who 
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often  provided  the  refreshing  distraction  that  helped  clear  my  mind  for 
the  task  at  hand.  I  would  most  like  to  express  thanks  to  my  Lord  and 
Saviour,  Jesus  Christ,  who  by  his  grace  has  faithfully  seen  me  through 
this  learning  experience. 

Alan  L.  Jones 
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Abstract 

The  purpose  of  this  3tudy  was  to  Investigate  the  suitability  of 
algorithms  derived  from  the  study  of  fractal  geometry  to  the  specific 
problem  of  image  segmentation.  The  use  of  two  widely  used  methods  and 
an  original  hybrid  technique  for  calculating  fractal  dimension  is 
demonstrated.  All  three  techniques  appeared  to  be  effective  to  some 
extent  in  segmenting  images,  but  shared  a  common  problem  with  the 
establishment  of  suitable  thresholds  for  robust  segmentation. 


vii 


IMAGE  SEGMENTATION  VIA  FRACTAL  DIMENSION 


I.  Introduction 


Overvlstt 

This  chapter  is  intended  to  provide  a  brief  discussion  of  the 
general  pattern  recognition  problem,  and  discuss  the  fractal  approach 
to  a  select  area  of  that  problem.  Encompassed  in  the  approach  is  a 
brief  introduction  to  the  fractal  set,  current  fractal  research  a3 
applied  to  pattern  recognition,  and  the  underlying  assumptions  used  by 
the  author.  Next,  a  plan  of  attack  is  developed,  the  goals  and  scope 
of  the  present  research  are  examined,  and  evaluation  criteria 
considered*  Finally,  the  support  equipment  used  to  carry  out  the 
research  is  discussed. 

The  Pattern  Recognition  Problem 

Pattern  recognition  is  of  interest  to  the  Air  Force  for  automated 
data  entry  and  reading  machines,  automated  recognition  of  individuals 
for  security  systems,  smart  targeting  systems  for  munitions,  as  well  as 
other  applications.  It  is  widely  accepted  that  the  only  effective 
general  purpose  pattern  recognition  machines  at  the  present  - .  .e  are 
the  visual  systems  of  living  organisms.  One  of  the  greatest 
difficulties  to  date  has  been  to  obtain  an  adequate  mathematical 


description  of  the  targets  of  interest.  To  paraphrase  Kabri sky's 

summarization  of  the  inherent  difficulty  in  pattern  recognition: 

To  find  arbitrary  targets  in  arbitrary  backgrounds  it  would  be 
useful  to  characterize  the  targets  mathematically.  In  practice, 
much  of  the  elegant  mathematics  developed  for  communication  and 
radar  systems  (detection  and  estimation  theory)  falls  apart  in 
many  realistic  scene  analysis  problems.  The  "signal"  implied  by 
the  shape  of  the  target  is  unknown  mathematically,  as  is  that  of 
the  background  noise.  As  a  result,  pattern  recognition  techniques 
are  quite  often  composed  of  a  combination  of  ad  hoc  3teps  that 
sometimes  work  depending  on  the  particular  data  set  and  individual 
problem  (10). 

An  examination  of  the  difficulties  facing  any  machine  tasked  with 
recognizing  images  brings  the  present  lack  of  consistent  success  into 
sharper  focus.  Closely  tied  to  this  inability  (in  general)  to 
adequately  describe  a  target  mathematically  are  the  following  five 
complications: 

1.  Translation  -  the  target  may  not  be  in  the  center  of  the 
sensor's  field  of  view,  and  in  general  is  not. 

2.  Scale  -  the  target  may  not  be  the  same  size  as  the  model 
"image"  that  the  pattern  recognizer  is  looking  for. 

3.  Rotation  -  the  target  may  have  an  In-plane  rotation  that 
differs  from  the  stored  model. 

U»  Aspect  -  the  target  say  havs  an  out-of-plane  rotation  that 
differs  from  the  stored  model  (for  example  the  sensor  is  viewing 
the  rear  of  the  target,  but  the  side  view  is  the  one  stored). 

5.  Noise  O0). 

Various  ad  hoc  approaches,  some  more  successful  than  others,  have 
been  used  to  overcome  these  problems.  One  of  the  more  successful 
approaches  was  an  algorithm  developed  by  Horev  (8)  and  further  expanded 
by  Kobe!  and  Kartln  (:2>.  Their  technique  was  shown  in  some  cases  to 
be  effective  in  dealing  with  the  translation,  scale,  and  (to  a  limited 
extent)  rotation  of  input  images  within  cluttered  scenes. 
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The  Fractal  Approach 


The  thrust  of  this  research  effort  is  to  apply  the  branch  of 

mathematics  known  as  fractal  geometry  to  the  pattern  recognition 

problem.  Mandelbrot  has  defined  a  fractal  set  to  be  a  set  for  which 

the  Hausdorff  Besicovitch  dimension,  D,  3trictly  exceeds  the 

topological  dimension  (14:15).  That  is,  fractals  have  noninteger 

dimension.  Consider  the  following  statement  by  Pentland: 

The  world  that  surrounds  ua,  except  for  man-made 
environments,  is  typically  formed  of  complex,  rough,  and  jumbled 
surfaces  ....  If  we  are  to  develop  machines  competent  to  caal 
with  the  natural  world,  therefore,  we  need  a  representational 
framework  that  is  able  to  describe  such  shapes  succinctly. 
...Fractal  functions  appear  to  provide  such  a  model,  in  part, 
because  many  basic  physical  processes  produce  fractal  surfaces 
(and  thus  fractals  are  quite  common  in  nature),  but  perhaps  even 
more  importantly  because  fractals  look  like  natural  surfaces. 
...This  is  important  information  for  workers  in  computer  vision 
because  the  natural  appearance  of  fractals  is  strong  evidence  that 
they  capture  all  of  the  perceptually  relevant  shape  structure  of 
natural  surfaces  (20:661,662). 

The  observation  that  fractals  may  be  an  appropriate  measurement 
feature  space  for  some  classes  of  pattern  recognition  problems  was 
precisely  the  motivation  for  this  research  effort. 

decent  Fractal  Research.  Fractal  analysis  has  been  applied  to  a 
surprisingly  diverse  number  of  disciplines  including  landsbre  analysis 
C14;25-33,2c:s?-276),  coscoleglc  distributions  (14:64-96),  surface  and 
volume  analysis  (14:109-115;  22),  medical  imaging  (13).  generation  of 
computer  graphic.-  (5:424-435;  1;  *6;  19;  23;  24:28-331,  image 
compression  (2),  modelling  of  metallisation  growth  on  semiconductors 
OS;  24:12;  25),  the  analysis  of  fractures  in  metal  and  stone  ( 14:460, 


and  aspects  of  scene  analysis  (10;  20;  21),  to  specifically  point  out  a 


Pentland's  research  has  demonstrated  that  fractal  analysis  can  be 
extremely  useful  in  segmenting  an  image  (dividing  it  properly  into 
areas  of  interest),  and  also  exhibits  an  ability  to  determine  whether 
or  not  the  segments  are  appropriately  modelled  as  fractals  (i.e. 
natural  objects)  (20),  Robust  segmentation  can  be  a  useful  first  step 
in  the  pattern  recognition  process,  depending  on  the  pattern 
recognition  technique  used  (10).  Robust  in  this  sense  implies 
repeatability  that  is  meaningful. 

Assumptions.  In  virtually  all  disciplines  to  which  fractal 
analysis  has  been  applied,  the  fractal  set  used  has  been  that  known  as 
fractional  Brownian  motion,  or  fBm.  This  research  effort  will  make  the 
assumption  that  the  fractal  set  of  interest  is  fBm.  This  fBm  set  ■ 
possesses  interesting  properties  that  particularly  suit  it  to  pattern 
recognition  in  natural  scenes.  These  properties  include  the  following: 

1.  A  three  dimensional  surface  with  a  spatially  isotropic  fractal 
Brownian  shape  produces  an  image  whose  intensity  surface  is 
fractal  Brownian  and  whose  fractal  dimension  is  identical  to  that 
of  the  components  of  the  surface  normal,  given  a  Lambertian 
surface  reflectance  function  and  constant  illumination  and  albedo, 

2.  A  linear  transformation  of  a  fractal  Brownian  function  is  a 
fractal  Brownian  function  with  the  same  fractal  dimension. 

3.  The  fractal  dimension  of  a  fractal  Brownian  function  is 
invariant  over  transformations  of  scale.  ( 20 : 66M— 665 ) 

The  first  of  these  properties  assures  that  any  fractal  analysis 

performed  on  the  two  dimensional  representation  (using  present  day 

input  sensors)  of  the  scene  stored  in  the  computer  hardware  is  valid, 

since  none  of  the  fractal  dimension  information  is  lost  by  acquiring  an 

image  rather  than  working  with  the  object  itself.  It  also  indicates 

this  form  of  scene  analysis  may  solve  the  pattern  recognition  aspect 

problem.  The  second  of  these  properties  predicts  that  fractal  analysis 


may  be  useful  in  dealing  with  the  problems  of  translation  and  rotation, 
since  both  are  linear  operations.  The  final  property  obviously 
indicates  that  fractal  analysis  will  yield  useful  pattern  recognition 
results  independent  of  scale.  Pentland  has  in  fact  demonstrated  this 
final  property  over  an  eight  to  one  scale  range  (20:668). 

Plan  of  Attack 

There  was  no  previous  AFIT  research  in  this  area,  so  all  of  the 
initial  groundwork  had  to  be  laid.  This  included  an  investigation  of 
fractals  in  general,  algorithms  to  determine  fractal  dimension, 
(Station  techniques  that  used  the  fractal  dimension,  and  the 
generation  of  software  to  accomplish  these  fundamental  tasks.  The  next 
step  was  to  apply  these  algorithms  to  natural  images  and  investigate 
the  performance  of  the  segmenter.  The  final  effort  was  to  apply  the 
refined  algorithm  to  realistic  Air  Force  images  to  investigate  its 
performance . 

Goals  and  Scope 

The  original  goal  of  this  research  was  to  find  a  fractal  kernel 
which  would  describe  segmented  regions  of  an  image  of  interest.  The 
end  result  would  have  been  a  "Fractal  Transform",  analogous  to  other 
familiar  transforms  (such  as  the  Fourier  transform),  possibly  rule- 
based  in  operation.  This  goal,  however,  proved  to  be  unrealistic.  The 
goals  which  were  finally  established  can  be  divided  into  three  broad 
areas.  The  first  was  to  generate  a  appropriate  fractal  decomposition 
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of  an  initial  set  of  images  which  would  be  used  as  standards.  This 
yielded  algorithms  potentially  capable  of  segmentation.  The  second 
goal  was  to  actually  segment  natural  images  in  a  robust  manner.  The 
final  goal  was  to  examine  typical  images  of  Air  Force  interest  to 
Investigate  the  performance  of  the  algorithm  in  realistic  image 
environments.  The  scope  of  the  research  was  limited  to  eight  images 
consisting  of  three  standard  images,  three  images  of  natural  scenes 
with  trees,  mountains,  etc.,  and  two  images  containing  possible 
military  targets.  Each  of  these  images  was  acquired  using  the  video 
digitizer  described  below.  Appendix  A  contains  each  of  the  images 
used. 

Evaluation  Criteria 

The  reference  for  evaluating  success  or  failure  in  this  effort  was 
the  human  visual  system,  specifically  the  author's  visual  system.  The 
achievement  of  each  of  the  above  goals,  and  each  algorithm's  relative 
merits,  was  determined  strictly  on  the  basis  of  the  following 
questions:  "Does  this  specific  technique  appear  (visually)  to  work?", 
and,  "How  does  It  compare  to  the  other  techniques  (visually)?" 

Support  Equipment 

The  primary  equipment  required  for  this  effort  was  a  MicroVAX  II 
AI  workstation  manufactured  by  Digital  Equipment  Corporation  (DEC). 

This  workstation  was  running  under  NicroVHS  and  was  equipped  with 

nine  megabytes  of  main  memory,  two  ?1  megabyte  hard  disk  drives,  the 
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GPX  color  upgrade  from  DEC,  a  FORTRAN  compiler  (VAX  FORTRAN  4.0),  and  a 
LISP  interpreter/compiler  (VAX  LISP  2.0).  There  was  also  an  occasional 
need  for  the  use  of  tha  video  digitizer  that  was  supported  by  another 
DEC  MicroVAX  II,  which  was  also  running  under  MicroVMS  4.4.  Both 
machines  were  linked  for  communication  purposes  via  the  DECNET  at  AFIT. 

The  digitizer  board  was  a  FG-100  series  singla  board  image 
processor  manufactured  by  Imaging  Technology  Incorporated  and  was  used 
to  digitize  images  with  eight  bits  of  gray  value  resolution  (256  gray 
levels)  per  picture  element  (pixel).  The  monochrome  images  were 
raptured  and  stored  as  files  by  the  ITEX  100  software  library  of  image 
processing  subroutines  supplied  by  Imaging  Technology  Incorporated.  A 
DADE  model  650  video  camera  equipped  with  a  Canon  6x18,  1 8-1 08mm,  f.^,5 
TV  zo«m  lens  was  us  id  as  the  source  of  video  input  to  the  digitize.'. 

Conclusion 

Pattern  recognition  is  an  inherently  difficult  task  for  a  machine 
to  perform  ac  the  present  time,  largely  due  to  the  inability  to 
succinctly  model  real  'jorld  images  in  a  mathematically  traotable  form. 
Compounding  this  are  tne  problems  of  translation,  rotation,  scale, 
aspect  changes,  and  noise  whose  solution  must  be  incorporated  into  the 
recognition  algorithm.  Fractals  appear  to  be  particularly  well  suited 
to  modelling  natural  scenes.  The  applicability  and  potential  of  this 
form  of  analysis  to  pattern  recognition  has  alrtidy  been  demonstrated. 
The  thrust  of  the  author's  research  was  to  first  verify  portions  of 
previous  fractal  analysis,  and  second,  to  extend  the  application  of 
this  analvsio  to  real  world  scenes  of  Interest  to  the  Air  Force. 
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II.  Fractals 


Overview 

This  chapter  is  intended  to  provide  the  necessary  background  in 
fractal  geometry  to  grasp  the  fundamentals  of  the  reasoning  behind  the 
algorithms  that  were  used  in  this  research.  A  general  introduction  to 
fractals  precedes  a  discussion  of  the  two  broadly  classed  types, 
deterministic  fractals  and  random  fractals.  The  various  concepts 
embodied  in  random  fractals  are  explored  with  an  eye  to  those  which  are 
particularly  relevant  to  the  algorithms  discussed  in  chapter  three. 

An  Introduction  to  Fractals 

It  has  only  been  recently  that  researchers  have  come  to  realize 
that  fractal  geometry  is  a  useful  tool  for  modelling  the  natural  world 
around  us.  Yet,  as  Mandelbrot  has  pointed  out,  fractals  do  in  fact 
suocinotly  model  many  natural  phenomenon  (14).  These  inolude,  but  are 
by  no  means  limited  to,  coastlines  (14:25-33),  galactic  clusters  (14:84- 
96),  craters  (14:301-309),  turbulence  (14:97-105),  flesh  (14;147-150), 
linguistics  (14:341-348),  trees  (14:151-165),  meteorology  (14:461), 
earthquakes  (14:461),  river  discharges  (14:247-255),  and  economics 
(14:335-340).  The  fractal  dimension  of  a  surface  is  also  an  excellent 
indicator  of  its  perceived  roughness  (20:662;  21:254). 
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The  standard  connotations  associated  with  dimension  can  be 
attributed  to  Euclid  (circa  300  B.C.)  who  began  his  Book  I  (plane 
geometry)  and  Book  XI  (spatial  geometry)  with  some  of  the  following 
fundamental  definitions: 

1.  A  point  is  that  which  has  no  part. 

2.  A  line  is  breadthless  length. 

3.  The  extremities  of  a  line  are  points. 

A  surface  is  that  which  has  length  and  breadth  only. 

5.  The  extremities  of  a  surface  are  lines. 

6.  A  solid  is  that  which  has  length,  breadth,  and  depth. 

7.  An  extremity  of  a  solid  is  a  surface.  (14:409) 

A  more  general  view  of  dimension  was  formulated  by  Hausdorff  and 
Besicovitch.  Hausdorff *s  contribution  was  the  use  of  test  functions 

j 

(14:364)  of  the  form  h(x)  =  g(d)x  to  place  a  mathematical  measure  on 
some  object  S.  For  instance,  consider  the  area  of  a  circle  and  its 
radial  measure.  In  this  case,  the  object  S  would  be  the  circle,  h(x) 
would  be  the  area,  g(d)  a  constant  (pi),  x  the  radius,  and  d  equal  to 
two.  For  Euclidean  objects,  the  "d"  in  any  Hausdorff  test  function 
will  be  an  integer  equal  to  the  object's  topological  dimension. 
Besicovitch  extended  this  concept  (14:364)  to  the  case  where  d  is  not 
an  integer,  and  S  is  not  a  standard  shape.  Besicovitch  went  on  to  show 
that  "for  every  set  S  there  exists  a  real  value  0  (corresponding  to  the 
d  in  Hausdorff's  test  function  equation)  such  that  the  d-oeasure  is 
Infinite  for  d<D  and  vanishes  for  d>D"  (14:364).  It  is  this  D  that  is 
the  fraotal  dimension  of  the  object  in  question. 

As  briefly  discussed  in  chapter  one,  a  fractal's  dimension  exceeds 
the  topological  dimension.  But  what  does  "fractal"  mean? 

I  coined  fractal  from  the  Latin  adjective  fractus.  The 

corresponding  Latin  verb  f range re  means  "to  break:"  to  create 

irregular  fragments.  It  is  therefore  sensible  —  and  how 

appropriate  for  our  needs!  —  that,  in  addition  to  "fragmented" 


(a3  in  fraction  or  refraction)  fractus  should  also  mean 
"irregular,"  both  meanings  being  preserved  in  fragment. 

The  proper  pronunciation  is  frac*  tal,  the  stress  being 
placed  as  in  frac1  tion.  (14:4) 

So  fractals  are  irregular,  fragmented,  and  usually  of  noninteger 
dimension.  In  fact,  sets  of  integer  dimension  may  or  may  not  be 
fractal,  but  every  set  that  has  noninteger  dimension  is  assuredly  a 
fractal  set  (14:15). 

A  better  understanding  of  this  fractional  dimension  can  be  had  by 
examining  the  concept  of  self-3imiliarity.  "When  each  piece  of  a  shape 
is  geometrically  similar  to  the  whole,  both  the  shape  and  the  cascade 
that  generate  it  are  called  self-similar"  (14:34).  Self-similarity  is 
the  foundation  upon  which  all  fractals  are  built,  whether  in  a  strict 
sense  as  is  the  case  for  deterministic  fractals,  or  in  a  statistical 
sense  in  the  case  of  random  fractals.  This  property  of  self-similarity 
exists  at  any  3cale  (theoretically),  and  for  this  reason  fractals  are 
often  referred  to  as  being  scaling  (24:4).  This  scaling  property  does 
not  hold  over  all  values  in  practical  applications,  but  does  hold  over 
a  wide  enough  range  to  be  extremely  useful.  This  range  varies  with  the 
application;  in  computer  generated  graphics  it  can  be  as  high  as  10^  to 

Q 

10  (24),  whereas  in  image  data  processing  it  may  only  be  as  high  as 
four  to  ten  (20:668).  Typically  the  range  is  constrained  by  the  field 
of  view  of  the  input  sensor,  the  resolution  of  the  input  sensor,  or  the 
resolution  of  the  output  device,  although  other  factors  may  contribute. 

Another  closely  related  concept  is  that  of  similarity  dimension. 

As  some  mathematicians  might  interpret  it,  the  similarity  dimension  and 
the  Hausdorff  Besicovitch  dimension  may  differ,  but  in  the  instances 
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where  they  agree,  it  is  also  the  fractal  dimension  (14:37).  This 
similarity  dimension  can  be  calculated  from  the  equation  Nr°  =  1 ,  or 
equivalently: 

D  =  log  N  /  log  (1  /r)  (1 ) 


where 


N  =  the  number  of  self-sis?  ..-ar  parts 
r  =  the  ratio  of  the  part  to  the  whole 
D  a  the  similarity  dimension 

Consider  how  this  applies  to  Figure  1.  In  the  case  of  the  self-similar 
line  segment,  N  is  five,  1/r  is  five,  and  D  is  log(5)  /  log(5),  or 
obviously  one.  In  the  case  of  the  self-similar  triadic  Koch  curve, 
attributed  to  von  Koch  1904  (14:35),  N  is  4,  1/r  is  3,  and  D  is 
log(4)  /  log(3),  or  approximately  1.2618.  This  is  in  fact  an  instance 
where  the  similarity  dimension  is  the  fractal  (Hausdorff  Besicovitch) 
dimension.  The  following  observation  by  Voss  may  provide  additional 
insight: 


As  D  increases  from  one  toward  two  the  resulting  "curves"  progress 
from  being  "line-like"  to  "filling"  much  of  the  plane.  Indeed, 
the  limit  (as)  D  (approaches]  two  gives  a  Peano  or  "space-filling" 
curve.  The  fractal  dimension  D,  thus,  provides  a  quantitative 
measure  of  wlggllness  of  the  curves.  Although  these  von  Koch 
curves  have  fractal  dimensions  between  one  and  two,  they  remain  a 
"curve"  with  a  topological  dimension  of  one.  The  removal  of  a 
single  point  outs  the  curve  in  two  pieces.  (24:6) 


Deterministic  Fractals 

With  minor  exceptions,  the  bulk  of  the  groundwork  necessary  for 
understanding  deterministic  fractals  was  established  in  the 
introduction.  "The  mathematical  principle  behind  the  generation  of 
these  (deterministic)  fractal  shapes  involves  the  iteration  of 
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Figure  1.  Standard  vs  Fractal  Self-Similarity  (14:44) 
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algebraic  transformations"  (16:62),  specifically  of  mapping  Euclidean 
space  onto  itself.  This  concept  of  mapping  is  the  underlying  principle 
uaed  in  the  computer  generation  of  deterministic  fractals.  As  might  be 
assumed,  recursive  programming  is  often  the  method  used  to  implement 
these  mappings  (23)*  To  this  point,  the  only  mapping  considered  was  in 
the  real  plane,  but  other  classes  of  deterministic  fractals  exist  that 
are  generated  by  mappings  in  the  complex  plane.  These  include  what  are 
known  as  Julia  sets  and  the  Handelbrot  set  (14:180-192;  19;  23).  with 
transformations  typically  of  the  fora  f(z)  a  z  -  u  where  u  is  a 
complex  number  and  z  may  or  may  not  be.  Figure  2  is  an  illustration  of 
an  example  from  the  Mandelbrot  set,  which  uses  a  mapping  in  the  complex 
plane.  Note  that  the  fractal  boundary  of  this  figure  is  quite 
different  from  the  triadie  Kook  curve  of  Figure  1*  These 
transformations,  whether  in  the  real  or  complex  plane,  are  applied  as 
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shown  in  Eq  (2)  (16:62).  Note  how  this  equation  clearly  illustrates 
the  inherently  recursive  nature  of  deterministic  fractal  generation. 

xQ  -  x,  =  T(x),  x2  =  TET(x)3,  x3  =  1 £TCT(x) ] } ,  ...  (2) 

where 

x_  s  the  initial  point  in  the  plane  (real  or  complex) 

Tu  =  the  functional  form  of  the  transformation 


Random  Fractals 

The  primary  difference  in  deterministic  and  random  fractals  i3  in 
the  interpretation  of  self-similarity.  Deterministic  fractals  look 
exactly  the  same  at  all  scales  ;  random  fractals  do  not.  Each  smaller 
Dortion  will  look  like,  but  not  exactly  like,  the  larger  object  from 
which  it  came  (24:7).  The  fractal  dimension  of  random  fractals  can 
still  be  calculated  from  Eq  (1)  and  formulations  similar  to  it,  but  now 
the  D  calculated  is  an  average  (24:7).  To  add  another  element  of 
complication,  random  fractals  may  also  exhibit  statistical  self- 
affinity,  rather  than  the  simpler  case  of  statistical  self-similarity. 
Under  an  affine  transformation,  the  mapping  described  by  Eq  (2)  can  now 
have  different  scaling  factors  for  each  of  the  Euclidean  coordinates 
(2:17).  Thus,  the  X  coordinates  may  be  scaled  differently  than  the  V 
coordinates,  which  in  turn  may  be  scaled  differently  than  the  Z 
coordinates,  for  example.  Compare  this  then  to  the  simpler  case  of 
statistical  self-similarity,  where  each  of  the  coordinates  would  be 
scaled  equally.  The  fractals  formed  by  such  an  affine  transforeatior 
are  properly  known  as  self-affine,  rather  than  self-similar,  but  the 
distinction  is  not  always  clearly  made  in  the  literature.  The  term 
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self-similar  is  often  used  in  a  generic  sense  that  incorporates  self¬ 
affinity  (His  349). 

As  pointed  out  in  chapter  one,  the  fractal  set  of  interest  in  tnis 
research  is  that  known  as  fractional  tor  fractal)  Brownian  motion 
(fBm),  and  it  is  this  subject  which  will  occupy  the  remaining  material 
in  this  chapter.  Strictly  speaking,  fBm  are  statistically  self-affine, 
which  can  cause  algorithms  for  calculating  fractal  dimension  to  produce 
ambiguous  results  depending  on  the  measurement  technique  used 
(24:19,20).  The  reason  is  simple:  the  scaling  may  vary  with  each 
Euclidean  dimension,  and  it  is  inevitably  these  scalings  which  are  used 
in  some  fashion  to  calculate  the  fractal  dimension. 

Fractional  Brownian  motion  has  its  basis  in  ordinary  Brownian 
motion,  first  described  by  R.  Brown  around  1827  while  observing 
particles  in  a  liquid  suspension  under  a  microscope  (23).  This  erratic 
movement  can  be  modelled  by  a  random  process  whose  increments  are 
Gaussian  distributed  (14:351;  23;  24:16;).  Additionally,  the  mean 
square  of  these  increments  has  a  variance  proportional  to  the 
difference  in  the  independent  variable  (23;  23416).  In  the  case  where 
the  independent  variable  is  time  (t),  this  can  be  more  succinctly 
stated  by  Eq  (3): 

<IX(t^)  -  X(t^ ) l2>  e  K  lt2  -  t, I  (3) 

where 

<  >  denotes  statistical  expectation 

K  denotes  a  proportionality  constant 
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Figure  3  is  a  typical  sample  of  one  dimensional  Brownian  motion. 

Eq  (3)  can  be  generalized  to  (24:16) 

<lX(t2)  -  X( t, )  I2>  =  K  lt2  -  t^2*1  (4) 

where 

<  >  denotes  statistical  expectation 

K  a  a  proportionality  constant 

H  a  the  Hur3t  exponent,  in  interval  [0,1]  (14:249) 

Eq  (4)  is  a  mathematical  generalization  of  Brownian  motion  in  one 

dimension  that  will  model  fBm.  When  H  is  equal  to  1/2,  the  formulation 

is  simply  that  of  ordinary  Brownian  motion.  The  cases  of  greatest 

interest  are  for  H  in  the  Intervals  [0,1/2)  and  (1/2,1).  In  the  former 

case,  the  increments  of  X(t)  are  positively  correlated;  and  in  the 

latter  case,  negatively  correlated  (24: 16).  The  Hurst  exponent  is 

directly  related  to  the  fractal  dimension  of  the  process  being  modelled 

by  the  relation  (24:24) 

D  *  E  ♦  1  -  H  (5) 

where 

D  a  the  fractal  dimension 
E  a  the  Euclidean  dimension 
H  =  the  Hurst  exponent 

The  effect  of  II  on  a  typical  trace  can  be  seen  in  Figure  4.  3y 
examining  Figure  4  and  knowing  that  fractal  dimension  increases  with 
perceived  roughness,  it  should  be  obvious  to  the  reader  that  a  larger  H 
corresponds  to  a  scalier  fractal  dimension,  which  is  precisely  what 
Eq  (3)  indicates. 
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Figure  3.  One  Dimensional  Brownian  Motion  (23) 
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Figure  5i.  The  Effect  of  H  on  f3a  (23) 


Conclusion 


When  it  comes  to  fractals,  there  really  is  "no  new  thing  under  the 
sun.  Is  there  any  thing  whereof  it  may  be  3aid,  See,  this  is  new?  it 
hath  been  already  of  old  time,  which  was  before  us"  (4).  As  has 
been  pointed  out,  objects  which  can  be  modelled  well  by  fractal 
geometry  pervade  the  world  around  us.  An  understanding  of  thi3 
geometry,  which  is  new,  moves  it  from  a  past  of  captivating 
mathematical  oddities  to  the  realm  of  a  truly  useful  analytical  tool  in 
the  present. 
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TII»  Experimental  Method 


Overview 

In  the  discussion  that  follows,  two  techniques  commonly  seen  in 
the  literature  for  determining  the  fractal  dimension,  D,  will  be 
reviewed,  as  well  as  a  original  hybrid  technique  investigated  by  the 
author.  The  two  types  of  threshold  processing  employed  and  the 
reasoning  behind  them  will  also  be  discussed.  The  actual  application 
of  these  algorithms  is  discussed  within  each  sub-section.  Appendix  B 
contains  the  actual  programs  used. 


Dimension  Measurement 

A  wide  variety  of  techniques  for  measuring  D  are  discussed  in  the 

literature.  Two  of  these  are  used  by  the  author,  and  are  apparently 

x 

the  most  commonly  used;  they  are  the  methods  employing  the  1/f  nature 
of  the  power  spectral  density  (PSD)  of  the  Fourier  transform  of  the 
signal  of  interest,  and  what  is  known  as  box  dimension.  Both  of  these 
techniques  are  particularly  well  suited  to  implementation  in  FORTRAN 
(25),  largely  due  to  its  array  handling  capabilities  and  the  highly 
optimizing  nature  of  the  compiler.  The  final  original  technique  used 
does  not  calculate  D  specifically,  but  does  borrow  from  concepts  used 
in  the  computer  generation  of  fBm.  It  was  investigated  simply  because 
of  its  ease  of  Implementation  and  relatively  fast  computation  time. 
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In  each  case,  the  general  approach  used  was  to  divide  the 
512  X  480  pixel  image  into  8X8  pixel  regions,  which  will  be  referred 
to  as  sub-images  throughout  the  rest  of  this  discussion,  for  subsequent 
processing.  This  size  of  sub-image  yielded  a  segmentation  resolution 
of  64  X  60  regions.  The  number  of  regions  could  have  been  increased, 
of  course,  by  simply  choosing  a  smaller  size  sub-image,  but  the  other 
consideration  was  the  Gaussian  nature  of  the  mean  square  increments  of 
fBm.  It  was  the  author’s  conviction  that  64  pixel  samples  should  be 
sufficient  to  insure  (within  a  reasonable  margin  of  error)  that  the 
increments  between  intensity  values  were  in  fact  Gaussian  in  nature, 
but  this  was  never  investigated.  The  three  fractal  dimension 
measurement  techniques  that  follow  were  applied  to  each  of  the  8X.8 
pixel  sub-images. 

PSD  Rolloff.  The  PSD  of  fBm  is  proportional  to  1 /f*  (24:24), 
where  x  can  be  related  to  the  Hurst  exponent  H  by  (24:24;  23): 

x  =  2H  ♦  1  (6) 

(For  a  development  of  this  relation,  see  Appendix  C.)  Substituting  for 
H  in  Eq  (5),  yields  a  relationship  between  the  spectral  rolloff  x,  the 
Euclidean  dimension  E,  and  the  fractal  dimension  D: 

D  :  E  ♦  (3  *  x)  /  2  (7) 

This  technique  was  implemented  by  making  use  of  the  Wlener- 
Khintohine  relation  (6:209),  which  simply  states  that  the  magnitude 
squared  of  the  Fourier  transform  of  a  signal  (i.e.  -  the  PSD)  is  a 
Fourier  transform  pair  with  the  autocorrelation  of  the  original 
signal.  It  was  felt  that  a  simple  X  and  Y  autocorrelation  would  be 
advantageous  from  the  viewpoint  of  computation  time  versus  implementing 
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an  8  X  8  two  dimensional  fast  Fourier  transform  from  which  the  PSD 
could  then  be  calculated.  Thus,  the  PSD  was  determined  using 
autocorrelations  in  the  X  and  Y  directions  (with  no  offset  in  Y  or  X, 
respectively)  for  each  of  the  8X8  sub-images,  which  were  then 
transformed  with  a  discrete  cosine  transform  to  obtain  the  spectral 
coefficients.  The  3pectral  rolloffs  in  X  and  Y  were  then  calculated  by 
simply  arithmetically  averaging  the  slope  in  the  seven  harmonics 
present  (the  DC  term  was  ignored),  and  these  two  were  then 
arithmetically  averaged  to  yield  a  rolloff  that  was  characteristic  of 
the  entire  8X8  pixel  sub-image.  This  was  then  used  to  calculate  the 
fractal  dimension,  D,  via  Eq  (7).  The  Euclidean  dimension  E  was 
assumed  to  be  three. 

Box  Dimension.  Consider  an  image  comprised  of  a  set  of  m  points, 
S,  in  E^dimensional  Euclidean  space,  which  a  priori  can  be  considered 
to  be  equiprobable,  and  with  associated  probablility  measure  P(ra,L). 
"P(m,L)  is  the  probability  that  there  are  m  points  within  an  E-cube  of. 
size  L  centered  about  an  arbitrary  point  [ml  in  S"  (Voss:2U).  This 
probability  measure  P(m,L)  is  then  normalized  such  that  the  summation 
over  all  the  points  m  in  S  is  equal  to  one  for  each  value  of  L 
(2*1:25).  The  mass  dimension  (or  box  dimension)  is  then  simply  the 
first  moment  with  respect  to  ra  of  P(m,L)  (2*1:25).  This  is  used  to 
calculate  the  fractal  dimension  D  from  the  relation  (2*1:25) 

M(L)  =  KLD  (8) 

where 

M(L)  =  the  mass  dimension  for  a  particular  value  of  box  size  l, 

K  =  a  constant  of  proportionality 

D  =  the  fractal  dimension 
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This  technique  was  implemented  by  first  normalizing  all  pixel  gray 
values  in  the  8X8  pixel  sub-image  to  the  interval  [1,83.  Thi"  can  be 
thought  of  as  yielding  a  cube  7  units  on  a  side  in  the  image  intensity 
space.  This  also  functioned  as  a  form  of  energy  normalization.  P(m,L) 
wa3  initialized  by  considering  it  as  a  discrete  distribution  of  64 
bins,  each  of  which  represented  the  number  of  points  within  the  cube, 
for  each  value  of  L.  There  were  64  bins,  because  for  large  enough 
values  of  L,  the  cube  would  contain  all  of  the  image  intensity  points. 
The  bins  were  then  filled  by  considering  each  point  in  the  sub-image  to 
be  the  center  of  a  cube  ranging  from  an  L  of  two  units  on  a  side  up  to 
an  L  of  seven  units  on  a  side,  and  counting  the  number  of  points  within 
that  cube.  This  yielded  the  six  distributions  P(m,2),  P(m,3),  ...» 
P(m,7)»  eaoh  of  which  was  then  normalized  to  one,  so  that  it  could  be 
treated  as  a  discrete  probability  density  function.  The  first  moment 
with  respect  to  m  was  then  calculated  for  each  P(m,L)  yielding  M(2), 
M(3),  ....  M(7).  A  value  proportional  to  the  fractal  dimension  D  of 
the  sub- image  was  then  calculated  using  Eq  (8)  by  calculating  the 
average  slope  of  a  log-log  mapping  of  L  vs.  M(L) . 

Hybrid  Brown.  This  technique  has  no  precedent  in  the  literature, 
although  it  borrows  from  concepts  used  in  the  computer  generation  of 
fBm  for  terrain  simulations  via  one  dimensional  Brownian  functions  (see 
reference  (23)).  In  this  method,  the  8x8  pixel  sub-image  intensity 
surface  is  normalized  in  X,  Y,  and  Z  to  [0,1).  Note  that  this  mapping 
is  substantially  different  from  the  past  techniques  in  which  the  X  and 
Y  coordinates  retained  their  mapping  on  the  Interval  [1,83,  and  their 
values  only  changed  in  integer  intervals.  The  sub- image  array  is  then 
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sorted  based  on  the  intensity  value  (the  Z  coordinate)  from  the  lowest 
value  to  the  highest  while  simultaneously  storing  the  path  taken 
through  the  array.  This  path  is  then  treated  as  three-space  Brownian 
motion,  and  the  increments  from  point  to  point  in  X,  Y,  and  Z  as  the 
path  is  traced  are  found.  This  yields  63  increments  in  each  of  the 
coordinates  X,  Y,  and  Z.  For  each  of  the  coordinates,  the  mean  square 
of  these  increments  is  then  calculated.  Since  the  assumption  is  that 
the  path  is  indeed  Brownian,  each  coordinate's  mean  square  value  can  be 
assumed  to  be  independent  from  the  other  two.  The  mean  square  that 
characterizes  the  8X8  pixel  sub-image  simply  becomes  the  addition  of 
the  three  previously  calculated  values.  The  Hurst  exponent  is  then 
calculated  via  Eq  (4),  where  the  independent  variable  interval  is 
assumed  to  be  1/63.  The  interval  chosen  might  appear  to  be  somewhat 
arbitrary,  but  there  are  63  increments  that  are  an  artifact  of  the  64 
element  sub- image,  and  it  seemed  fitting  that  an  interval  of  1/63  was 
appropriate.  It  is  important  to  note  that  what  is  being  calculated  in 
this  technique  is  not  the  fractal  dimension  of  the  sub- image;  it  is 
merely  a  statistical  calculator  based  on  fractal  Brownian  motion.  It 
was  investigated  because  of  the  relatively  quick  computation  time 
involved,  and  to  compare  its  performance  in  its  ability  to  segment  to 
the  known  fractal  dimension  calculators. 

Threshold  Techniques 

The  natural  problem  that  arises  once  the  calculators  finish  is  how 
to  organize  the  data  into  some  graphically  meaningful  output.  Two 
approaches  were  tested.  In  each  case,  the  fractal-based  measurement 
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was  used  to  determine  which  color  (or  equivalently,  intensity)  bin  a 
sub-image  would  be  mapped  to.  The  sub-images  were  mapped  to  four 
different  values. 

Linear.  This  method  was  implemented  by  scaling  the  working  array, 
which  was  the  6*1  X  60  array  containing  the  output  of  the  calculators, 
to  range  in  value  from  zero  to  four.  Thresholds  were  then  set  at  one, 
two,  and  three,  and  each  of  the  corresponding  elements  in  the  array  was 
then  assigned  a  value  based  on  which  of  the  thresholds  the  region  was 
between . 

Adaptive.  This  method  used  histogram  manipulation  to  determine 
where  the  minima  were  in  a  256  bin  histogram  of  sub-image  values,  which 
in  turn  were  used  to  set  the  array  thresholds.  It  was  implemented  by 
first  scaling  the  working  array  to  values  in  the  interval  [0,2561.  The 
number  of  values  in  the  working  array  in  the  interval  (255,2561  was 
then  placed  into  bin  255,  the  number  of  values  in  the  interval 
(25^,2551  into  bin  25^,  etc.  The  histogram  was  then  numerically 
integrated  by  simply  replacing  each  bin  value  by  the  running  total  sum, 
and  an  average  slope  for  the  entire  integrated  histogram  was  then 
calculated.  The  slope  from  bin  to  bin  was  then  compared  to  the 
average.  Any  slope  that  was  greater  than  the  average  was  considered  to 
be  due  to  either  a  peak  or  the  rising  edge  of  a  peak  in  the  original 
histogram.  This  approach  allowed  the  computer  to  search  for  maxima, 
and  indirectly,  the  minima  as  well.  The  search  proceeded  from  the 
largest  intensity  values  downward  until  three  minima  were  found.  These 
bin  values  between  peaks  then  corresponded  to  the  three  thresholds  that 
were  set  for  assigning  working  array  elements  to  their  segmented 
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value.  These  threshold  values  were  simply  an  artifact  of  the  natural 
clusters  of  fractal  dimension.  The  hope  was,  of  course,  that  the 
clusters  corresponded  to  what  the  human  eye  would  have  segmented  in  the 
original  image. 

Modified  Linear.  This  method  can  be  thought  of  as  scaling  the 
working  array  into  eight  levels.  The  lower  five  levels  were  all  mapped 
to  black,  the  next  higher  level  to  a  dark  gray,  the  next  higher  level 
to  a  light  gray,  and  the  upper  level  to  white.  This  is  still  a  linear 
method  for  establishing  thresholds,  it  simply  ignores  the  lower  four 
values.  This  technique  emphasises  variations  in  the  upper  half  of  the 
values  in  the  working  array  histogram,  and  essentially  ignores  those 
variations  in  the  lower  half. 

Conclusion 

The  implemenation  of  a  fractal  dimension  calculator  is  not 
intuitively  obvious.  Two  standard  techniques  were  presented  along  with 
a  discussion  of  the  specific  implementation  used  by  the  author. 
Additionally,  an  original  technique  and  the  reasoning  behind  it  was 
discussed.  Finally,  the  three  methods  used  for  establishing  thresholds 
and  the  implementation  of  each  was  discussed. 


IV.  Results 


Overview 

This  chapter  examines  the  success  of  the  three  techniques  for 
calculating  fractal  dimension  as  applied  to  the  segmentation  of  the 
eight  images  mentioned  in  chapter  one.  Each  technique  is  discussed 
separately  with  accompanying  examples.  Additional  results  can  be  found 
in  Appendix  D. 


Output  Comments 


All  three  methods  of  establishing  thresholds  were  found  to  be 
inadequate.  The  adaptive  threshold  never  successfully  set  the 
thresholds  as  it  was  hoped.  This  was  perhaps  due  to  the  Jaggedness  of 
the  histogram  profiles.  The  linear  method  was  in  almost  every  case  an 
absolute  failure.  The  most  successful  of  the  thresholding  methods  was 
the  modified  linear  method.  After  examining  histograms  of  the 
distribution  of  fractal  dimension  within  various  scenes,  the  reason 
became  obvious.  It  appears  that  the  fractal  dimension  calculators  do 
in  fact  differentiate  between  different  types  of  natural  objects,  as 
evidenced  by  the  clusters  observed  in  the  histograms.  The  problem  with 
this  is  that  the  clusters  vary  In  number,  location,  and  with  the 
fractal  dimension  calculator  used.  This  was  precisely  the  reason  that 
the  linear  method  failed  completely.  Additionally,  these  clusters 
tended  to  be  within  the  upper  half  of  the  histogram  which  explains  why 


the  modified  linear  threshold  technique  met  with  some  measure  of 
success.  All  graphic  results  were  obtained  using  the  modified  linear 
threshold  technique. 

PSD  flolloff 

This  technique  appears  to  work  from  histogram  distributions  and  is 
relatively  fast,  typically  processing  an  image  in  about  two  minutes. 
This  technique  does  not  separate  the  clusters  within  the  histogram 
enough  to  yield  visually  pleasing  results  with  the  present  threshold 
approaches.  This  can  be  readily  seen  by  considering  Figure  6,  which  is 
the  result  of  the  simple  input  image  of  Figure  5,  as  opposed  to  the 
result  of  Figure  7,  which  is  the  result  of  the  more  complex  image  of 
Figure  8.  Note  that  intuitively  Figure  5  would  contain  one  prominent 
cluster  in  the  histogram,  but  Figure  7  would  contain  several. 

Box  Dimension 

With  the  present  threshold  methods,  this  technique  seems  to  be  the 
most  effective  in  segmenting  an  image,  but  is  very  slow,  typically 
processing  an  image  in  *i5  minutes  or  more.  Figure  9  is  an  example  of 
the  effectiveness  of  this  approach  on  the  image  shown  in  Figure  8. 

This  teohnique  tends  to  spread  out  the  locations  of  the  clusters  within 
the  histogram,  and  thus  is  more  likely  to  separate  them  with  the  fixed 
thresholds  established  by  the  modified  linear  method  of  thresholding. 
This  teohnique  was  also  U3ed  without  first  scaling  the  intensity 
values,  yielding  virtually  identical  results. 


■■j-.g-i.  sima  it 


Figure  5.  Input  Inage  (14:68) 
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Figure  6.  Segaentatloe  via  PSD  Rolloff 


© 


4-tJ 


Hybrid  Brown 


This  technique  also  appears  to  work,  although  like  the  PSD 
technique,  it  does  not  separate  the  clusters  within  the  fractal 
dimension  histogram  as  well  as  the  box  dimension  technique.  It  is  also 
relatively  fast,  typically  processing  an  image  in  about  two  minutes. 
Figures  10  and  11  are  examples  of  the  output  of  this  technique  with 
Figures  5  and  8  as  input  images,  respectively.  Note  the  similarity  of 
these  results  to  that  of  the  PSD  rolloff  technique. 


V,  Conclusions  and  Recommendations 

Conclusions 

Some  general  conclusions  can  be  drawn  from  the  results  presented 
in  the  previous  chapter.  One  of  the  most  important  is  that  the  use  of 
modelling  natural  objects  using  fractal  geometry  shows  great  promise  in 
image  processing.  The  use  of  fractal  geometry  to  mathematically  model 
regions  of  interest  in  an  image  is  a  fundamentally  different  approach 
than  the  Euclidean  models  and  Euclidean-extracted  measurement  spaces 
used  in  pattern  recognition  to  date.  Of  the  three  techniques 
investigated,  all  three  appear  to  work  to  some  extent.  The  PSD  rolloff 
technique  and  hybrid  Brown  technique  both  fail  in  cluttered  images, 
apparently  due  to  an  inability  at  the  present  to  adaptively  set  the 
necessary  thresholds  to  provide  robust  segmentation.  All  of  the 
techniques  appear  capable  of  extracting  useful  image  information  from 
noisy  (real  world)  images,  the  box  dimension  technique  demonstrably  so. 

Recommendations 

The  following  general  areas  of  investigation  appear  to  be  worth 
pursuing.  The  foremost  area  of  investigation  is  that  of  adaptively 
setting  the  thresholds,  such  that  each  of  the  techniques  examined  in 
this  researoh  would  be  useful  in  terms  of  robust  segmentation.  Two 
possible  approaches  to  this  might  be  searching  for  minima  in  a  curve 
that  has  been  3pline-fit  to  the  histogram  profile;  or  secondly, 


correlating  a  generic  peak  with  the  histogram  and  thresholding  the 
result  to  determine  the  location  of  the  maxima  (and  thus  the  minima). 
Secondly,  investigate  the  distribution  of  intensity  increments  within 
each  of  the  processed  sub-images.  If  in  fact  regions  have 
distributions  that  are  not  Gaussian,  theory  indicates  the  region  is  a 
likely  candidate  to  be  part  of  a  man-made  object,  simply  because  the 
model  for  fBm  does  not  fit  the  data.  More  simply  stated,  natural 
objects  fit  the  fBm  model,  man-made  objects  do  not. 

Areas  of  investigation  concerning  the  specific  implementation  of 
the  .software  as  it  presently  exists  are  those  of  pre-processing  and 
region  resolution.  Preprocessing  the  image  may  aid  the  dimension 
calculators  in  achieving  the  desired  end.  No  edge  enhancement,  or  any 
other  form  of  pre-processing,  was  employed  before  passing  the  image  to 
the  fractal  dimension  calculators.  Also,  is  it  in  fact  necessary  to 
use  8X8  pixel  sub-images?  There  is  probably  some  optimum  size  of  sub¬ 
image  that  is  a  good  tradeoff  between  processing  speed  and  processed 
image  resolution. 

Turning  to  the  original  goal  of  this  research,  which  was  to  find  a 
fractal  transform,  one  area  of  current  research  has  been  encountered 
whloh  may  further  progress  on  this  problem.  The  general  approach  would 
be  based  on  an  extension  of  the  work  by  Barnsley  discussed  in  reference 
(2).  Barnsley’s  efforts  focus  on  finding  a  set  of  graphic  kernels 
(that  can  be  modelled  by  fractal  geometry)  which  can  be  used  to 
reconstruct  an  entire  object.  The  general  approach  is  to  use  a  fixed 
set  of  probabilistic  affine  transformations  known  as  iterated  function 
set  codes  to  generate  this  image.  The  image  can  be  thought  of  as 
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closely  representing  the  object  in  a  mean  square  error  sense.  If  this 
approach  could  be  automated  (it  is  not  presently)  and  applied  to  an 
input  image  that  has  been  partitioned  on  -  grid,  it  may  be  possible  to 
find  a  fractal  kernel  that  represents  each  partition.  The  fractal 
dimension  of  this  kernel  can  then  be  compared  to  the  fractal  dimension 
of  the  entire  partition  using  the  calculators  developed  in  this 
research  in  an  iterative  fashion  until  some  predetermined  measure  of 
fit  is  obtained.  This  approach  may  or  may  not  work,  contingent  on 
whether  it  is  a  true  assumption  that  the  fractal  dimension  of  the 
kernel  used  to  generate  a  partition  is  equal  (or  proportional)  to  the 
fractal  dimension  of  the  partition  itself.  An  issue  not  discussed  by 
Barnsley  is  the  orthogonality  of  such  fractal  kernels,  or  of  any 
effects  that  might  be  caused  by  discontinuities  in  the  fractal  kernel 
from  partition  to  partition.  It  is  not  intuitively  clear  that  any  such 
decomposition  would  yield  orthogonal  kernels  similar  to  the  harmonics 
used  in  a  Fourier  transform.  It  is  also  quite  possible  that  there  are 
in  fact  discontinuity  phenomenon  similar  to  the  Gibb’s  phenomenon 
exhibited  by  discontinuities  in  a  signal  which  is  Fourier  transformed. 

The  final  recommendation  is  specifically  targeted  to  anyone 
pursv-ing  future  work  with  fractal  geometry.  The  author  considers  it  a 
necessity  that  at  a  minimum  references  (5),  (20),  (23),  and  (2H)  be 
3'.udled  and  thoroughly  understood.  These  references  will  provide  a 
solid  foundation  for  further  study  in  the  application  of  fractal 
geometry  to  image  processing. 
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Appendix  A:  Input  Images 

This  appendix  simply  displays  in  numerical  order  the  images  that 
were  used  in  this  investigation.  The  file  in  which  the  image  was 
stored  is  given  for  each.  These  images  were  acquired  with  the 
digitizer  and  video  camera  described  in  chapter  one.  The  files  were 
converted  into  a  format  readable  by  the  FORTRAN  oompiler  using  the 
ITEX_TO_RMS  program  written  by  Capt.  Richard  Roberts. 
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ototyp*  McDonnell  Douglas  F-1SC  fitted  with  conformsl  fuel  tanks  anti  with  a  LANTIRN  navigational 
pod  undar  its  port  angina  air  intaka 


I  external  stores  stations  remain  available  with  the  CRTs 
use.  and  McDonnell  Douglas  has  developed  for  the  F-15 
tew  weapon  attachment  system  which  can  extend  the 
erating  radius  with  >arge  external  loads  by  up  to  40  per 
it.  Known  as  tangential  carriage,  it  involves  the  instal- 
ton  of  rows  of  stub  pylons  on  the  lower  comer  and 


MS1P  improvements  include  a  tactical  electronic  system 
consisting  of  a  Northrop  Enhanced  ALQ-135  internal 
countermeasures  system.  Lora!  ALR-56C  radar  warning 
receiver.  Tracor  ALE-4S  chaff  disperser,  and  Magnavox 
electronic  warfare  warning  system.  A  further  $274-4  million 
contract  was  received  in  December  1983.  Flight  testing  of 
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Appendix  B:  Software 


This  appendix  is  3imply  a  listing  of  the  various  routines  that 
were  used  to  investigate  image  processing  using  fractal-based 
algorithms.  Included  are  the  programs  that  ./ere  used  to  convert  the 
original  output  of  the  ITEX  digitizer  board  (the  only  C  program),  load 
the  image  files  into  the  VAX  FORTRAN  environment,  display  images  on  the 
GPX  hardware,  and  all  of  the  processing  routines.  Many  thanks  to  Capt. 
Richard  Roberts  and  Capt.  Dennis  Ruck  for  their  help  in  adapting  the 
ITEX  conversion  program  to  suit  the  author's  needs. 
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This  subroutine  calculates  the  rolloff  in  the  PSD  of  the  FFT 
spectrum  by  calculating  the  autocorrelation  in  X  and  Y* 


SUBROUTINE  AUTO ( ARRAY , SUB  ARRAY , X_INDEX , Y  INDEX) 

*  Define  variables  ... 

* 

INTEGERM  X  INDEX,  Y  INDEX 
REAL**  SUB  ARRAY,  CORR  X,  CORR  Y 
DIMENSION  CORR  X(0:7),~CORR  Y(0:7) 

DIMENSION  ARRAY(64 ,60)  ,  SUB  ARRAY (*, 8) 

*  — 

*  Loop  over  each  shift  ... 

* 

DO  1010  K-0,7 

4 

*  Zero  the  correlation  vectors  ... 

* 

CORR  X ( K )  -  0.0 
CORR~Y(K)  ■  0.0 

* 

*  Loop  ovir  «v%ry  of  th#  aub-i»»g#  . . . 

4 

DO  1000  J-1,8 

DO  1300  I#K+1 , 8 

4 

*  Autocorrelation  across  X  ... 

4 

CORR_X(K)  -  CORR  X(K)  +  SUB  ARRAY  ( I ,  J  ) 

1  *  SUB_ARRAY  ( I~X ,  0  ) 

1000  CONTINUE 

DO  1010  J-K-tl.e 
DO  1010  1*1,8 

« 

*  Autocorrelation  across  Y  ... 

4 

CORR  Y ( K )  -  CORR  Y ( K )  *  SUB  ARRAY (I , J ) 

1  *~SU8  ARRAY (iTJ-K) 

1010  CONTINUE 

4 

*  Pass  the  correlation  values  to  the  Fourier  routine  ... 

4 

CALL  DCT(CORR_X) 

CALL  PCT ( CORR  Y) 

DO  1020  K*0,7" 

4 

*  Add  small  value  to  keep  log  function  from  crashing  ... 

4 

CORR  X ( K )  *  CORR  K ( K )  ♦  0.000031 
CORR~Y ( X )  «  CO«R“Y<K)  ♦  0.000001 

*  Take  the  log  of  each  harmonic  ... 

4 

CORR  X ( X )  ■  DLOG10(CORR  X  Ik ) ) 

CORR~YtK)  »  OLOOl0tCORR_Y<K) ) 

1020  CONTINUE 


Calculate  the  rolloff  in  the  PSD  in  log-log  fora  ... 

Initialise  alope  values  to  aero  ... 

SLOPS  X  »  0.0 
SLOPEJf  -  0.0 

Calculate  slope  in  each  increment  and  add  togsthar  ... 
DO  1030  1*1.6 

SLOPE  X  »  SLOPE  X  ♦  (CORR  XII)  -  CORR  X  t  X  *  1 ) >  / 

1  ( ALOGlO ( FLOAT ( 1 1  I  -  A LOG  10  I  FLOAT (1*1) I  I 
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SLOPE  Y  ■  SLOPE  Y  +  (CORK  Y { I )  -  CORR  Y(I+1))  / 
1  (ALOOIO  (  FLOAT) X ) )  “  ALOOiO ( FLOAT ( 1  +  1 ))  ) 

1030  CONTINUE 


Crudely  calculate  th«  slop*  of  th*  lln*  by  averaging  .?. 

SLOPE  X  -  SLOPE  X  /  6.0 
SLOPE Jf  -  SLOPE~Y  /  6.0 

Calculate  th*  fractal  dinension  in  X  and  Y  ... 

DIM  X  -  3.0  +  0.5  •  (3.0  +  SLOPE  X) 

DIH~Y  -  3.0  +  0.5  •  (3.0  SLOPE~Y) 

Average  and  put  rasult  in  th*  work  array  ... 

ARRAY (X  INDEX, Y  INDEX)  -  (DIM  X  ♦  DIM  Y)  /  2.0 

RETURN 

END 


This  subroutine  siaply  swaps  two  valu*s  (us*d  for  sorting) 


SUBROUTINE  SWAP(X,Y) 

D*fin*  variables  ... 

REAL* 8  X,  Y,  TEMP 

Parfora  th*  swap  ... 

TEMP  -  X 
X  «  V 
Y  *  TEMP 
RETURN 
END 


This  subroutine  calculates  diaensionality  by  treating  th* 
inage  pixel  values  as  Brownian  notion  ss  th*  'particle* 
tree**  out  a  three  spec*  path  over  ( 1 0 , 1 ) , ( 0 , 1 )  ,  ( 0 , 1 1 )  ... 


SUBROUTINE  BROWN ( ARRAY , SUB_ARRAY , X_I NDEX , Y_I NDEX ) 

Define  variable*  ... 

LOGICAL* 2  ERROR 
INTEGER**  X  INDEX, Y  INDEX 
REAL**  ARRAY 

REAL* 8  SUB  ARRAY.  SORT,  DELTA,  SUM,  SQUARE 
DIMENSION  ARRAYtS*. 60) ,  SUB  ARRAY (8,8),  DE LTA (63) 
DIMENSION  SQUARE ( 3  I ,  SORT(ST.I),  SUHtS) 

Nornelie*  th*  array  SUB_ARRAY  to  nln  *  0,0,  *•»  »  1.0  ... 

CALL  HORN (SUB_ARRAY. ERROR) 

Check  error  trap  on  return  Iron  NORN  ... 

IP  (ERROR  .CQ.  .TRUE.)  THEN 

ARRAY) X  INDEX, Y  INDEX!  -3.0 
RETURN 
END  IP 


Now  associate  n'raslited  X,  Y,  end  S  coordinate*  so 
that  th*  eetu*:  ^*th  will  b*  known  ... 

DO  1000  J-1,8 
DO  1900  1*1,8 

SORT ( (J-J)*l»I,i)  .  tl-l!  /  7.0 
SORT)  (O'*  11*1*1,2!  -  )J-ll  /  7.0 


B-3 


SORT ( ( J-l ) *8+1 • 3 )  -  SUB  ARRAY (X,J) 

1000  CONTINUE 

Sort  ascending  2  values  , . . 

DO  1010  1-1,63 

Initialise  2  ... 

VALUE  -  SORT (I, 3) 

Test  each  of  tba  2  values  ... 

DO  1010  J-I+1,64 

ir(SORT(J,3)  .LT.  VALUE)  THEN 
VALUE  -  SORT (J, 3) 

CALL  SWAP ( SORT (1,1) , SORT ( J , 1 ) ) 

CALL  SWAP ( SORT (1,2), SORT (0,2)) 

CALL  SWAP ( SORT (1,3), SORT ( J , 3 ) ) 

END  If 

1010  CONTINUE 

*  Loop  through  X,  Y,  2  ... 

OO  1020  J-1,3 
SUM ( J )  -  0.0 
SQUARE ( J )  -  0.0 

* 

*  .  Oet  the  increaants  ... 

* 

DO  1020  1-1,63 

DELTA!  I)  -  SORT  (  X  +  l  ,  J  )  -  SORT(I.J) 

*  Calculate  eua  and  square  values  for  the  variance  ... 

« 

SUM ( J  )  -  SUM ( J )  ♦  DELTA ( I ) 

SQUARE ( J )  -  SQUARE ( J )  ♦  DELTA(X)  •*  2.0 
1020  CONTINUE 

*  Calculate  the  variance  in  X,  V,  and  2  ... 

« 

VAR_X  -  (63.0  •  SQUARE ( 1 )  -  (SUM(l)  ••  2.0))  /  3906.0 

VAR~ Y  -  (63.0  •  SQUARE ( 3 )  -  (SUM(2)  ••  2.0))  /  3906.0 

VAR- 2  »  (63.0  •  SQUARE ( 3 )  -  ( SUM ( 3 )  •»  2.0))  /  3906.0 

*  Aeeuaing  Brownian  notion  iaplies  that  the  variances  in 

*  X,  V,  and  S  can  be  treated  as  independent  Oeussian  r.v.'s 

VAR  -  VARJE  ♦  VAR_Y  ♦  VAR_l 

*  Now  calculate  H  (arbitrarily  actuate  delta  T  is  1/63)  ... 
H  -  ALOGl 0 ( VAR )  /  (2,0  '  ALOOIO t 6 3 . 0 ) ) 

*  And  place  in  array  ... 

ARRAY (X  INDEX .  Y  INOEXI  «  3.0  -  ABS(H) 

RETURN 

END 


This  subroutine  converts  the  work  array  into  an  image  stray. 


SUBROUTINE  CONVERT ( l N . OUT ) 

BYTE  OUT 
REAL'  4  !** 

DIMENSION  IN ( 1 4 , 60  ) ,  OUTIS12.440  ) 

Piece  Into  6X4  blocks  in  the  array  OUT  ... 

OO  1090  J-V.tO 
DO  1900  1-1.44 


H-4 


n 


DO  1000  K-1,8 

DO  1000  L»l,8 


1000 


1000 


continue 

RETURN 

END 


OUT ((I-l)*8+L,( J-l )  *  8  +  K )  ( IN ( 1 ,  J )  -  128) 

» 


1010 


1020 


This  subroutine  parforas  a  discrete  cosin*  transfora  on  ths 
autocorrelation  to  b*  used  in  calculating  tha  rolloff  of  th* 
power  spactral  dansity. 


SUBROUTINE  DCTIIN) 

REAL'S  IN,  PSD 
DIMENSION  IN ( 0 : 7 ) ,  PSD(7) 

Initialica  PSD  . . . 

DO  1000  1-1,7 
PSD(I)  -  0.0 
-CONTINUE 

Calculat*  tha  Pouriar  transforia  ... 

DO  1010  K-1,7 
DO  1010  N-0,7 

PSD(K)  -  PSD ( K )  a  IN ( N )  *  DCOSD(D8LE(N*K)M5.0) 
CONTINUE 

Load  raault  back  into  tha  array  to  ba  passad  ... 

DO  1020  1-1,7 

IN ( I )  -  FSD(I) 

CONTINUE 

RETURN 

END 


This  subroutin*  iaagai  an  array  of  bytaa  in  a  writing 
nod#  that  uaaa  th*  byta  valua  to  determine  tha  color 
(or  shad#  of  gray)  of  th*  pixal . 


SUBROUTINE  DRAW(IN  ARRAY  , COLOR , HEADER  ) 

BYTE  IM  ARRAY 

LOGICAL”!  COLOR 

CHARACTER* 12  HEADER 

REALM  RED,  GREEN,  BLUE,  INTENSITY 

DIMENSION  IH_ARRAY{512,4I0) 

DIMENSION  RED(  256I ,  GREEN  (  236 )  ,  BLUE ( 2  58 ) ,  INTENSITY t  256  ) 

Include  th*  UIS  graphic*  routine*  ... 

INCLUDE  *  SYS $ LIBRARY  tUISUSRDEP ' 

Craata  a  virtual  color  nap  ... 

VCN^IO  -  UIS$CREATE_COLOR_MAP(25« 1 

Croat*  a  virtual  display  and  associate  the  eoloraap  ... 

VO  to-  UIS5CREATE  DISPLAY (0. 0,0. 0,102*. 0,8*4. 0.16-1, IS. J. VCH_IO) 
CALL  UlSSGNA»LE_DlIPLAY_LISTt VOID) 

Check  to  toe  at  color  or  ehada*  of  gray  ... 

IP  (COLOR  .60.  .PALSE.I  THEN 
DO  1000  X-0,235 

•lack  to  whit*  in  ISC  atop*  ... 

INTENSITY ( I )  -  PLOAT(X)  /  255.0 
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1000 


CONTINUE 


*  rill  the  color  nap  ... 

» 

CALL  UXS$SET  INTENSITIES (VD  ID , 0 , 256 , INTENSITY ) 

ELSE 

*  Initialise  all  color  vactors  to  contain  0.0  ... 

* 

DO  1010  I  -  1,256 
RED ( I )  -  0.0 
GREEN ( I )  -  0.0 
BLUE! I)  -  0.0 
1010  CONTINUE 

* 

*  rill  thaa  with  color  ... 

* 

DO  1020  X-1,65 

VAL  -  0.30  +  ( FLOAT ( 1 )  *  (0.69  /  #5.0)) 

BLUE ( I )  *  VAL 
GREEN) #5  ♦  I)  -  VAL 
RED( 170  +  I)  -  VAL 
1020  CONTINUE 

« 

*  Don't  forget  white  ... 

* 

RED ( 256 )  -  1.0 
GREEN ( 256 )  -  1.0 
BLUE ( 256 )  -  1.0 

* 

*  Fill  the  color  nap  ... 

• 

CALL  UXS$SET  COLORS (VD_ID, 0,2 55, RED, GREEN, BLUE) 

END  ir 

CALL  UIS$SET_WRITINO_HODE ( VD_ID , 0,1, UIS$C_MODE_COPY ) 

*  Display  the  input  iaaga  array 

CALL  UXS$ IHAGE ( VD_ID , 1, 0. 0, 0.0, 512.0, 4#0. 0,512, 480, 8, IN_ARRAY ) 

*  Window  the  iaaga  ... 

« 

WO_ID  -  UX5SCREATE_WINDOW( VD_ID. * S YSSWORK3TATION ' .HEADER) 

*  This  subroutine  works  with  a  easll  bug  ... 

*  ...  to  gat  the  header  and  buried  windows  eove  the 

*  pointer  to  the  aide  of  the  iaaga  and  push  the  left 

*  aouse  t  :tt  3i\ 

*  select  EXIT  froa  the  aenu  ... 

RETURN 

END 


This  subroutine  is  used  to  read  in  the  ITKX  digitiser  files. 
The  files  Must  have  been  previously  converted  by  the 
I TEX  TORHS  executable  file  to  an  RMS  forest  (identified 
by  tKe  Tssi  extension  in  the  directory  ( IMAGES | )  to  be 
used  by  load. 


SUBROUTINE  LOAD  I  ARRAY , INTI LC ) 

BYTE  ARRAY 
CHARACTER* 32  INFILE 
DIMENSION  ARRAY ( 512.110) 

TYPE  1000 

1000  FORMAT!//'  lasge  f i V eXAMt . EXT?  VI 

*  Read  in  the  desired  lasge  tile  nsae  ... 

ACCEPT  1010,  INFILC 
1010  FORMAT  t  A3  2 ) 

iNriLC  -  ’ ihs v/infile 


*  IM  is  •  logical  to  access  tha  image  filas  on  tha  SMV2A  ... 

*  Connact  iaage  fila  to  unit  1  ... 

*  . 
OPEN ( 1 ,  FILE-INFILE,  READONLY , 

1  RECORDTYPE-'VARIABLE' , 

1  STATUS** 'OLD' ,  ACCESS* ' SEQUENTIAL' ,  FORM- ' UNFORMATTED '  ) 

* 

*  R«ad  in  th#  i»tg«  ... 

. 

DO  1020  J-l , 4 SO 

READ  (1)  (ARRAY(I.J) .1-1,512) 

1020  CONTINUE 
CLOSE! 1) 

RETURN 

END 


This  subroutine  calculates  dimensionality  using  what  is 
known  as  box  (or  mass)  dimension. 


SUBROUTINE  MASS ( ARRAY , SUB  ARRAY , X  INDEX, Y  INDEX) 

#  “  "* 

*  Define  variables  ... 

. 

LOGICAL* 2  ERROR 

INTEGER* 4  X  INDEX,  Y  INDEX,  LO_X,  LO  Y,  HX_X,  HI  V 
REAL* S  SUB_ARRAY,  MASSBIN,  M 

DIMENSION  MASSBIN! 7, 64 ) ,  H<7>,  ARRAY ( S 4 , 60 ) ,  SUB  ARRAY ( « , 8 ) 

*  Eero  all  of  the  elements  (bins)  in  MASSBIN  and  M  ... 

* 

DO  1000  1-1,7 
Mil)  -  0.0 
DO  1000  J-l, 64 

MASSBIN! I, J)  -  0.0 
1000  CONTINUE 
* 

*  Count  points  within  L  x  L  X  L  box  throughout  the  region  ... 

*  Pirtt  loop  incr»M«ntt  bon  slst  ... 

DO  1010  b-3,9 

*  t  two  loop*  incr«atnt  ovtr  ••ch  point  in  region  . . . 

DO  1010  o-v,« 

OO  1010  I  -  1.4 

*  Box  1 i.ite  in  t  ... 

* 

MIN  S  *  SUB  ARRAY (1,41  -  L 
KAX_S  -  iUBJARRAY! I ,J)  ♦  L 

*  »0#  lialte  in  V  ... 

• 

to  f  a  J  »  L 
sr“(t0  Y  . LT .  U  THEN 
W  V  •  I 
END  S? 

HI  ¥  «  J  *  L 
ir“(HI  ¥  .GT.  •)  THtN 
K)  T  •  I 
CRD  IF 

* 

•  Box  liaita  in  X  ... 

. 

LO  I  »  I  »  L 
I F~ I  LG  X  . LT .  1  I  THEM 
LO  I  -  1 
END  IF 
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V-iv,*, 


ft 

ft 


HI  X  -  I  +  L 
If“(HI  X  .GT.  8)  THEN 
HI  X  -  8 
END  I  r 

9 

Count  tha  points  in  tho  box  _ 

DO  1010  K-LO  Y , HI  Y 
DO  1010  N-LO  xThI  X 

IF  ( isub  arrayTn.k)  .LT.  HAX_E) 

1  .AND.  (SUB  ARRAY (N,K)  .GT.  MIN  Z ) )  THEN 

massbin(l-T, ; j-i ) *8+1 )  - 

1  MASSBIN( L-l , ( J-l ) *8+1 )  *  1.0 

END  IF 

Incraaant  L,  l,  a,  K,  N  ... 


1010  CONTINUE 

DO  1040  L-l, 7 


* 


1020 

1030 


1040 


1050 

10«Q 


Noraalisa  tho  boss  distribution  to  1.0  for  ooch  valuo  of  L  ... 

DO  1020  1-1,64 

M(L)  -  MIL)  +  MASSBIN ( L , I ) 

CONTINUE 
DO  1030  1-1,64 

MASSBIN ( L , X )  -  HASS BIN ( L, Z )  /  MIL) 

CONTINUE 

Calculata  tha  Bass  diaanaion  for  aach  L  ... 

DO  1040  1-1,64 

M(L)  -  MIL)  ♦  II  *  MASSBIN I L , I ) ) 

CONTINUE 

Now  aatiaata  tha  axpactad  valua  of  D  for  aach  L  ... 

DO  1050  L-l, 6 

MIL)  -  I  DLOG1 0  I  HI  L+l )  )  -  DLOG10  (MID)  )  / 

1  I DLOG10 I DBLC ( L+l ) )  -  DLOGIO I OBLEI L) ) 1 

CONTINUE 
D  -  0.0 
DO  1060  L-l, 6 
D  -  D  ♦  MIL) 

CONTINUE 

Plata  caleulatad  D  in  working  array  ... 


ARRAY  I X  INDEX. Y  INDEX  I  *  O  '  4.0 

RETURN 

END 


This  aubroutlna  alaply  noraalltas  tha  array  sub  to 

contain  8  v*luaa  batwaan  0  and  l.  aata  an  arrot  flap  if  all 

$  valuaa  ata  aqual. 


SUBROUTINE  NORM! SUB. ERR) 
REAL'S  SUB.  MM.  MX 
DIMENSION  IUBII.il 
LOGICAL*  2  ERR 


*  tnitiallaa  ain,  ■#«.  and  EBB  ... 

MN  -  25S.O 
MX  »  0.0 
ERR  «  -PAL2C. 


Rind  tha  Bin  and  aaa  E  valuaa  (ta  noraalita  and  aa 
an  arror  trap) 

DO  1000  J-l.S 
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ft**1 


r.‘ 


v*\*v  -v  r  Wif.y.'f,' 


.  XV.V.% '  V. 


1010 


DO  1000  1-1,8 

XP  ( SUB ( X , J )  .LT.  MN)  THEN 
HN  -  SUB { X , J ) 

END  IP 

XP  ( SUB ( X , J )  .OT.  HX)  THEN 
HX  -  SUB(I,J) 

END  IP 
CONTINUE 

This  is  ths  srcor  trap  ... 

IP( JXDXNT(MN)  .EQ.  JIDINT(HX))  THEN 

Slaply  sat  ths  srror  <1*9  ERR  »nd  axit  ... 

ERR  -  .TRUE. 

RETURN 
END  IP 

Noraaliaa  thw  array  SUB  to  -in  ■  0.0,  >ax  -  1.0 

RANGE  »  NX  -  MN 
DO  1010  J-1,8 
DO  1010  1-1 , 8 

Raaevs  tha  bias,  seals  values  ... 

SUB(I,J)  -  { SUB ( X , J )  -  HN)  /  RANGE 
CONTINUE 
RETURN 
END 


This  aubroutins  proesasss  ths  sagaantad  output  array  into  4 
bit  pattarns  (to  ssulats  ahadaa  of  gray)  which  ars  than  sant 
to  tha  DEC  LNO  J  plus  lasar  printac. 


SUBROUTINE  PRINT ( AERAY . THRESHOLD ) 

BYTE  DCS,  ST.  PF.  UP.  CR 
CHARACTER* J  CRLF 
CHARACTER* 5  IKIT 

CHARACTER’S  WHITE ,  GRAY*.  GRAY  1 1 ,  BLACK.  SPACE 
INTEGER’S  THRESHOLD 
REAL* 4  ARRAY 

DIMENSION  ARRAY ( 64 . (0  I .  THRESHOLD) J) 

Xnittallc#  prints?  control  byts*  and  graphics  strings  ... 


DCS  -  *  SO  *  X 
ST  »  ‘»C‘X 

tr  *>  ■  A'x 
CR  •  'O'K 

rr  -  *  c  ■  x 
CRtr  -  •$**• 
XMIT  -  *S(.’0gf 
WHITE  -  *im?) 
CRAY*  -  'fQJOfQ 
GRAY 1 1  «  •Tltm 
BLACK  - 
SPACE  -  • 


iGcaphict  »cd# 

(Graphic*  data  t-r-inator 
IPo#  t  inch  top  aargin 
(Standard  earriaga  taturfi 
(Porn  faad  to  ajact  ptga 
(Graphic*  carriage  rstuen  and  tici*  (sad 
(Initial!**  graphics  charactariat ica 
(Sixsl  graphics  string 


(Nova  iaaga  to  aiddta  af  ptga 


Fco-pt  t«  snaura  prints?  le  tsady  ... 

pause  ‘Prints?  switch**?  continue  <cr>.  or  exit  <ca>  ...  < 

Opan  tha  print-?  ... 

OPEN  < ACCESS-’ SEQUENT! At ‘ .  FORK-’ UNFORMATTED1 .  PtLS-‘£$A0: 
Cantar  lasgs  ta  ba  printad  ... 

WRITE  t  S )  LP. LP. UP. LP.LP.LP.CR, SPACE. SPACE 
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*  Plica  printar  in  graphics  aod«  ... 

WRITE! 3)  DCS,  INXT 

*  Begin  to  send  ectuel  graphics  data  ... 

* 

DO  1000  J-1,60 
WRITE! 3)  CRLP 
DO  1000  1-1,64 

IP  ( ARRAY ( I ,  J  )  . LE ,  THRESHOLD! 1 ) )  THEN 
WRITE! 3)  BLACK 

ELSE  IF  ( ARRAY ( I , J )  .L£.  THRESHOLD! 2 ) )  THEN 
WRITE! 3)  GRATIS 

ELSE  IP  { ARRAY ( I, J)  .LE.  THRESHOLD! 3 ) )  THEN 
WRITE! 3)  GRAY4 

ELSE 

WRITE! 3)  WHITE 
END  IP 

1000  CONTINUE 

*  Tamimta  graphics  and  ajact  paga  — 

WHITE* 3)  ST  p  ft 

*  Close  printer  and  exit 
« 

CLOSE! 3) 

RETURN 

END 


This  suhroutina  is  us.d  to  aanipulat  tha  work,  array  in 
a  visually  usefiui  torn. 


SUBROUTINE  PROCESS (IN, VALLEYS) 

LOG I CAL *2  ERROR.  CHANGE,  TEST 

CHARACTER*!  answer 

INTEGER*!  BIN.  VALLEYS 

REAL* 4  IN.  MN,  NX ,  RANGE.  SLOPE 

DIMENSION  IN!fi4.«QI,  BIN(0:25SJ,  VALLEYS ( 1 ) .  TBST144.601 

Inicialiaa  uln  and  a«c  ... 

HN  -  355.0 

NX  «  0.0 


*  Pind  tha  ain  and  was  valuaa  ... 

♦ 

DO  1000  J-i.SO 
DO  1000  I-5..44 

IF  (!!»(!. 3!  .  VT  •  NR)  .AND.  USm.J)  ,  RE .  O.OH  THEN 

«N  »  INU.J) 

END  IP 

IP  UNII.JJ  ,<5T.  NX)  VHfcB 
NX  -  XNtt.3) 

END  IP 

{900  CONTINUE 

*  WetNalita  tha  array  !H  te  win  »  0-0.  *»s  »  255.0  ear raapssdlng 

*  t*  the  gray  valuaa  .a  tha  GPX  . . . 

* 

RANGE  -  HS  -  It* 

DC  SOU  3*  1,49 
DO  1010  t- 1 .44 

*  Rataava  tha  biaa,  acala  valuaa  ... 

* 

IP  t  SN( S • J I  .HE.  0.0)  TWER 

SN.l.J)  •  2«<  0  *  |S«i:,3l  *  HR  f  /  SARDS 

END  IP 


n-io 


1010  CONTINUE 

VALLEYS ( 1 )  -  127 

VALLEYS! 2)  •  159 

VALLEYS ! 3 )  -  221 

CALL  S>RXNT(  IN,  VALLEYS) 

RETURN 

END 


This  subroutine  writes  t.'.e  work  array  to  disk. 


SUBROUTINE  WRITER ILE ( ARRAY , OUTPILE, TV PE ) 

Define  veriablee 

LOGICAL* 2  TYPE 
CHARACTER *1  CHARRAY 
CHARACTER»32  OUTPILE 
REAL* 4  ARRAY 

DIHEN.' ION  ARRAY  (64,60)  ,  CHARRAY  (44,6-0) 


Associate  OUTPILE  with  a  VHS  filespec  ... 
OUTPILE  *  ‘ DUA1 t { A JONES . DATA  1 '//OUTPI 1 E// ' -DAT* 


Create  an  outfits  and  associate  it  with  a  logical  unit  ... 
OPEN  t  PILE -OUTPILE .  PORK" *  FORMATTED' ,  STATUS- ' NEW  ,  UNIT  -  2) 
Write  the  work  array  to  disk  ... 


1090 

1910 


1020 


1010 

1040 


IP  (TYPE  .EQ.  .TRUE.)  THEN 
DO  1010  J-i.SO 

DO  1010  I-1.ST.4 

*  roue  of  4  reels  for  each  row  of  the  work  array  ... 

WRITE! 2.  1000)  (ARRAY  I $*N.2) .R-O.T) 
PQRKAYI 0P1Q , 4 ) 

continue 

ELSE 


Write  out  character  graphics  of  44  character*  for 
each  row  of  the  w-ock  array  ... 


DO  1020  2* 1.60 
Dvt  1020  l—l  ,  4.4 

IP  ( I  IK? 1  ARRAY! I .2) 1  CO. 

CHARRAY! I, 2)  «  •  1 

ELSE  IP  !  UKT!  ARRAY!  !  .  2  t  > 
CHARRAY! Z, J!  - 
ELSE  IP  (UNTt  ARRAY!  1.2)) 
CHARRAY (1,2!  -  *W 
ELSE 

CHARRAY! I. 21  »  **• 

END  IP 

CONTINUE 
00  10*0  2-1. 

WRITS (2,1515)  S CHARRAY S ! .2). 5*1 
PORMAT 1 4  4A1 I 
CONTINUE 
ENG  If 


0)  THEN 
.  C(J .  I)  THEN 
. E0 .  2 )  THEN 


.44) 


Dene  with  file  -»  »e  cleee  it  ... 

CLOSE! J ( 

RETURN 

EH.' 
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PROGRAM  HAIR 
INCLUDE  'GLOBAL. FOR' 

CALL  LOAD! IMAGE, STRING) 

TYPE  3 

rCRMATl/’  Frocaiaing  iaaga  ...  work  aeray  lino  •  '/) 

DO  10  J«l,60 
TYPE  5 , J 
FORMAT! 1*1 
DO  10  1*1 ,64 

Soqaont  iaaga  into  •  X  8  aub-iaagas  ... 

CALL  SUB  (IMAGE,  SU8_IHAGE.  I  *8-7,  .1*8-7) 

Calculcto  fractal  diaanaion  of  aaeh  via  hybrid  Brown  . 

CALL  BROWN (WORK, SUB  IMAGE, t,J) 

CONTINUE 
FLAG  *  .FALSE. 

Oiaplay  original  iaaga  ... 

CALL  DRAW! IMAGE , FLAG , STRING } 

PAUSE 

Proc»»»  fractal  diaahaioa  raauita  and  duap  to  printar 

CALL  PROCESS (WORK, THRESHOLDS) 

Convert  to  an  iaaga  file  (ia  -  to  bytaa)  ... 

CALL  CONVERT (WORK. IN AGE  5 

Taka  a  look  at  final  iaaga  ... 

CALL  ORAMUMaGE. FLAG. STRING) 

PAUSE 

END 


8-12 


y**************************************************************************** 


Program  name  ITEX_to  RMS.c  ? 

Author  Richard  Roberts 

Special  thanks  to  Stavva  Johnson  and  Rick  Rainas  for  assistanne* 


This  program  inputs  a  file  from  disk  into  an  array.  The  image  * 

storad  by  the  ITEX100  software  is  read  in  and  converted  into  * 

decimal  format.  Each  pixel  is  stored  into  the  two  dimensional  * 

array  with  a  color  value  between  0  and  255.  * 

* 

The  program  then  writes  the  image  to  dick  in  the  .SSI  format.  * 

The  program  assumes  images  of  si*e  512  by  480.  * 

ft 

. . . . . . . . . * . / 


^include  ”sys$l ibra ry : stdio .h" 
((include  "sys$library :  rms  .  h” 

int  i,j,k; 

char  name (2561: 

char  test (21;  /* 

char  char  headet(631; 

int  num_Keader l 6 3  1  ; 

char  comment [  256  ]  ; 


/»  loop  control  variables  */ 

/*  name  of  image  file  */ 
left  or  right  image  select  Vjt  */ 

/*  array  containing  header  char  info  */ 
/*  array  containing  header  num  info  */ 
/*  array  containing  user  image  comment 


*/ 


/*  change  these  array  values  to  match  ITER  image  */ 

char  im_char_array l 512 ! ( 480 1 ;  /*  array  containing  character  pixel  values  */ 

int  im_num_array ( 512 ] ( 480 1 ;  /*  array  containing  decimal  pixel  values  */ 

FILE  *afile,  * f open  (); 


((define  RECORD_SIZE  (siseof  record) 

int  rms_status; 
struct  FAB  rab; 
struct  RAB  rab; 


/•  change  row  data  array  if  not  using  512  pixel  rows  */ 

struct 

( 

char  row  data  ( 512  1 ; 

)  record; 


char 

•filename; 

int 

index ; 

int 

lest; 

cha  r 

response ; 

main!  ) 

( 

printf ( "\n\nCONVERT  ITEX100  Tormat  to  ,SSI\n\nM: 
response  *  ‘ V ' : 

while  (  ( response*- ‘ V ‘  )  |  ( response** ' y 1  )  ) 

( 

get  i n  f  o (  )  ; 
rea3_disk (  ) ; 

/*  Initialise  BHS  file  •/ 

/*  Change  file  extension  to  SSI  */ 

last  *  at rlen ( 4nam* )  -  1; 

name (last-21  »  'S': 

name ( l as t-  1  |  »  'S'; 

name ( 1 est  1  «  ' I ' ; 


fab  *  cc$rm»  fab; 


u-n 


*  «  v  A  ■  »  ** 


V  V  % 


f tb. f ab$l_fna  ■  knama; 
tr b .  f  ab$b_f ns  «  st rlan ( inaaa ) ; 
fab.  fab$b-org  -  FAB$C_SEQ; 
fab . f ab$b~rat  ■  FAB$M_CR; 
f ab . 2ab$v~»rs  «  0;  “ 

rab  ■  cc$r»s_rab; 
rab . rab$l_f ab  »  tfab; 

r»a_status  -  svs$craata  (afab); 

r«»_status  »  sys$connact  (irab); 

/*  changa  loop  countars  in  not  512  X  480  itaga  nota  "<"  in  loop  */ 

fortjaO;  j<480;  j++)  /*  writa  iaaga  pixal  colors  •/ 

{ 

for(i»0;  i<512;  i+a) 

{ 

racord . row_dat a t i 1  ”  in  nu»  array(i)(j)  -  128; 

) 

rab,rab$b  rac  -  RAB$C_SEQ; 
rab . rab$ l^rbf  «  iracord; 
rab .  rab$w_rar,  »  RECORD  SIZE; 
ras_statu¥  •  ayoSput  (Irab); 

} 

/*  Closa  RMS  fila  •/ 
ras_status  ■  syaSclosa  (ifab); 

prlntf ( "Convart  Anothar  Fila?  (Y/N)  '); 
scanf("%s“,trasponsa); 

};  /•  and  whila  •/ 


) 


/ . . . . . . 

gat  inf o ( ) 

(  " 

c  l  s  (  ) ; 

printf(“\n  Entar  tha  naaa  of  tha  fila  (includa  ,EXT)\n>“); 
scant ( "»a‘ . (naaa) ; 

raturn; 

1 


. . . . . . . 

raad  dial)  t  ) 

( 

afila  topanlnana,  *  r  ‘  ;  /*  opan  fila  for  raad  only  */ 

prlntf > *\n\n\n\n\n\n  Plaaaa  wait  whila  fila  is  procaaaad . \n  *); 

fortk-O;  k<«63;  kaa)  /*  raad  haadar  info  V 

( 

char  haadar(k)  «  gatctatilal; 
ou»_haada r ( h |  *  char  haadar(k); 
iftnua  haadar I k  |  <0  1 

nut  haadat|k)  «  nua  haadartkl  «  254:  ✓*  corract  tor  wrap  around  */ 

) 


if (nun  HaSdar ( 1 2 ((  ou»  haadar 1 1 2 |»»2 t  /•  chack  fila  foraat  typo  »/ 

I 

c  1 1  f  )  : 

prinlft*\n  Fila  can  not  ba  procass.d.*): 

prlntfl*\n  Fila  foraat  »u»t  ba  alghr 

pr int f ( '\n\n\n\n^n  Prat*  RETURN  to  guit.Nn.*!: 

gatcha  t  (  I 

gatchar  l  1  ■ 

raturn; 
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.MV, 


) 


£ o r ( ;  j <nua_heade r ( 2 ] ;  j++) 
consent! j ]~«  getc(afile); 


/*  reed  consent  irtt  */ 
» 


/*  change  loop  counters  if  not  512  X  480  iaago  —  note  "<"  in  loop  */ 

/*  read  iaage  pixel  colors  */ 


for(j»0:  j<480;  j++) 

( 

for(i-0;  i < 5 1 2 ;  i++) 

( 

ia_char_array I i 1 ( j )  -  getc(afile); 
i»”nun_array ( i ][ j )  *  i»_char_array t i )( j 1 ; 
ifTia_nua_array( i H  j  J  <07 

ia”nua  arraytiilj]  *  ia_nua_array ( i ] ( j ]  +  256; 
record. row  detail]  •  ia  nua  array(i](j]  -  128; 


) 


) 


f close (af 11* ) ; 
return; 


) 


/ . . 

cUI  ) 

l 

printf ( "\n\n" )  ; 
return; 

) 
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Appendix  C:  Fractal  Dimension  via  PSD  Rolloff 


This  appendix  is  simply  a  derivation  of  the  relationship  described 
by  Eq  (6).  It  is  taken  in  its  entirety  from  reference  (23). 

Consider  the  following  fundamental  property  of  fBm:  If  X(t) 
denotes  fBm  with  Hurst  exponent  0<H<1 ,  then  the  properly  scaled 
function 

Y(t)  =  r“HX(rT) 

for  a  given  r>0  has  the  same  statistical  properties,  and  thus  the  same 
spectral  der**lty,  as  X(t).  From  this  basic  property  and  the  use  of 
some  elementary  calculus,  the  results  of  Eqs  (6)  and  (7)  can  be  deduced 
as  follows. 

Let  us  fix  some  r>0,  and  let 

,  .  „  Y(t)  =  r”HX(rT) ,  0<t<T 

*  ”  0  ,  otherwise 

and  adopt  the  notation 

F„( f ,T) ,  Fy(f,T)  Fourier  transforms  of  X(t,T),  Y(t,T) 

3y ( f , T ) ,  sJ(f,T)  Spectral  densities  of  X(t,T),  Y(t,T) 

Sj|(f),  Sy(f)  Spectral  densities  of  X(t),  Y(t) 

We  then  compute  the  Fourier  transform  of  Y(t,T) 

FvCf,T)  =  =  r~H  f  X(s)e~2"i(r/r)3ds/r 

1  U  ■'O 

where  we  have  substituted  s/r  for  t  and  ds/r  for  dt  in  the  second 
integral.  It  follows  then,  that 

Fy(f,T)  =  r‘(H^1)Fx(f/r,rT) 

Thus,  for  the  spectral  density  of  Y(t,T) 

Sy(f,T)  =  r’‘(2H4',)(!/rT)lFx(f/r,rT)f2 
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-’i-a ■» ».  •’Lt »  k’.  y n sr*  V'T-*’?  y.  ^-v  y.- y*  ys 


and  in  the  limit  as  T  approaches  infinity,  or  equivalently  as  rT 
approaches  infinity,  we  conclude 

SY(f)  =  r"^2H+1)Sx(f/r) 


But,  since  Y  is  just  a  properly  scaled  version  of  X,  their  spectral 
densities  coincide  such  that 


Sx(f)  =  r'(2H+1)Sx(f/r) 


Now  we  formally  set  f=1  and  replace  1/r  again  by  f  to  obtain  the 
desired  result 


Appendix  D:  Output 


This  appendix  contains  the  graphic  output  using  the  modified 
linear  method  of  thresholding.  Each  figure  is  labelled  to  indicate  the 
fractal  dimension  calculator  and  input  image  (refer  to  Appendix  A) 
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